%This script implements the test ignoring estimation uncertainty that is reported in Table 1
%Corresponding author: Rodrigo Adao
%Date: 09/11/2024
%Input: model_data_fgkk.mat
%Output: est_noparest.mat 

%% Preliminaries
clear;
close all; 

%Define vector of large varieties in our sample
trade_target = 0.9;

%set lacal paths
function_path = "..\Functions\";
addpath(function_path,'-frozen');
set_local_paths;

%cluster: 1 - no cluster, 2 - sector-country 
test_cluster_grid  = 1;

Nspec = length(test_cluster_grid);
results_dW = zeros(13,3,Nspec);
results_dM = zeros(13,3,Nspec);
results_dT = zeros(13,3,Nspec);
results_dX = zeros(13,3,Nspec);

%% Estimation

for j=1:Nspec

tic
    
%% Import and adjust data

display('------------- running new j -----------------')
    j

%% Import data

display('Implement preliminary steps')

%estimation parameters
test_cluster = test_cluster_grid(j);

%Run preliminary steps for estimation
estimation_preliminary_steps;

clear    rho_* share_pm share_px share_rm delta_ig_0 a_star_ig_0 z_star_ig_0 a_ig_0 aM_g_0 AM_s_0 AD_s_0 betaT_0 beta_s_0 alphaL_s_0 alphaI_s_0 alpha_ks_0 barL_rs_0 Z_rs_0 ...
             m_ig_0 x_ig_0 tau_ig_0 tau_star_ig_0 Lsr Dsg m_val_sample m_val_adj adj_m_s adj_m_gs x_val_sample x_val_adj adj_x_s adj_x_gs

%% IV for welfare

 %Naive IV matrices
    var_adj=sum(ww_dW.^2)/Nobs;
    
    %No GE
    share_NC_nGE_dW  = shareIV_tau;    
    share_wNC_nGE_dW  = ww_dWmat_naive*share_NC_nGE_dW ;
    coef = MC_term1*share_wNC_nGE_dW;
    share_wMC_nGE_dW  = share_wNC_nGE_dW- C*coef;
    share_wMC_nGE_dW  = share_wMC_nGE_dW/var_adj;
    
%Auxiliary matrices for control and welfare adjustments
%load auxiliary adjustment vectors
load(save_share_path + "estimation_shares_t" + trade_target + ".mat", 'adjGE_wMC')

    adjGE_wMC_mat= spdiags(adjGE_wMC, 0 , Nobs, Nobs );
    aux = adjGE_wMC_mat*shareMOD_dW;
    share_wMC_dW_dW = aux - C*(MC_term1*aux);
    
    clearvars shareIV_tau aux adjGE_MC_mat coef ww_dWmat ww_dWmat_naive share_NC_nGE_dW share_wNC_nGE_dW adjGE_wMC_mat

%% Implement test

%Predictions 
    dx = shareMOD_dW*shifters;

    coef = MC_term1*dx;
    dx_MC_dW = dx - C*coef;
    dx_MC_dM = dx_MC_dW(indM==1);
    dx_MC_dT = dx_MC_dW(indT==1);
    dx_MC_dX = dx_MC_dW(indX==1);

    coef = MC_term1*dx(sample_naive);
    dxn_MC_dW = dx(sample_naive) - C*coef;
    dxn_MC_dM = dxn_MC_dW(indMn==1);
    dxn_MC_dT = dxn_MC_dW(indTn==1);
    dxn_MC_dX = dxn_MC_dW(indXn==1);    
    
%Outcomes
    dy = dWout_sample;

    coef = MC_term1*dy;
    dy_MC_dW = dy - C*coef;
    dy_MC_dM = dy_MC_dW(indM==1);
    dy_MC_dT = dy_MC_dW(indT==1);
    dy_MC_dX = dy_MC_dW(indX==1);

    coef = MC_term1*dy(sample_naive);
    dyn_MC_dW = dy(sample_naive) - C*coef;
    dyn_MC_dM = dyn_MC_dW(indMn==1);
    dyn_MC_dT = dyn_MC_dW(indTn==1);
    dyn_MC_dX = dyn_MC_dW(indXn==1);    
    
%Difference between outcome and predictions    
    dyt = dy - dx;

    coef = MC_term1*dyt;
    dyt_MC_dW = dyt - C*coef;
    dyt_MC_dM = dyt_MC_dW(indM==1);
    dyt_MC_dT = dyt_MC_dW(indT==1);
    dyt_MC_dX = dyt_MC_dW(indX==1);

    coef = MC_term1*dyt(sample_naive);
    dytn_MC_dW = dyt(sample_naive) - C*coef;
    dytn_MC_dM = dytn_MC_dW(indMn==1);
    dytn_MC_dT = dytn_MC_dW(indTn==1);
    dytn_MC_dX = dytn_MC_dW(indXn==1);    
    
clear C MC_term1 shareMOD_dW coef

%% Welfare stacked specification

disp('running estimation for welfare stacked outcomes')

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dW, dx_MC_dW]);
    cor_MC_dW = correl_pred(2,1);
    
%Test based on naive IV
    zn_wMC_dW = share_wMC_nGE_dW*shiftersIV;  
    
    [btt_zn_wMC_dW , set_zn_wMC_dW , rjt_zn_wMC_dW , r0t_zn_wMC_dW, pjt_zn_wMC_dW , p0t_zn_wMC_dW ] = implement_test_cluster(dytn_MC_dW, zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, critical, 0, cluster_vec);
    [bty_zn_wMC_dW , sey_zn_wMC_dW , rjy_zn_wMC_dW , r0y_zn_wMC_dW, pjy_zn_wMC_dW , p0y_zn_wMC_dW ] = implement_test_cluster(dyn_MC_dW , zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dW , sex_zn_wMC_dW , rjx_zn_wMC_dW , r0x_zn_wMC_dW, pjx_zn_wMC_dW , p0x_zn_wMC_dW ] = implement_test_cluster(dxn_MC_dW , zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    zw_wMC_dW = share_wMC_dW_dW*shiftersIV;  
    
    [btt_zw_wMC_dW , set_zw_wMC_dW , rjt_zw_wMC_dW , r0t_zw_wMC_dW, pjt_zw_wMC_dW , p0t_zw_wMC_dW ] = implement_test_cluster(dyt_MC_dW, zw_wMC_dW , share_wMC_dW_dW , shiftersIV, critical, 0, cluster_vec);
    [bty_zw_wMC_dW , sey_zw_wMC_dW , rjy_zw_wMC_dW , r0y_zw_wMC_dW, pjy_zw_wMC_dW , p0y_zw_wMC_dW ] = implement_test_cluster(dy_MC_dW , zw_wMC_dW , share_wMC_dW_dW , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dW , sex_zw_wMC_dW , rjx_zw_wMC_dW , r0x_zw_wMC_dW, pjx_zw_wMC_dW , p0x_zw_wMC_dW ] = implement_test_cluster(dx_MC_dW , zw_wMC_dW , share_wMC_dW_dW , shiftersIV, critical, 0, cluster_vec);

%Output
    
    results_zw_dW = full([bty_zw_wMC_dW, btx_zw_wMC_dW, btt_zw_wMC_dW; ...
                          sey_zw_wMC_dW, sex_zw_wMC_dW, set_zw_wMC_dW; ...
                          pjy_zw_wMC_dW, pjx_zw_wMC_dW, pjt_zw_wMC_dW; ...
                          p0y_zw_wMC_dW, p0x_zw_wMC_dW, p0t_zw_wMC_dW]);

    results_zn_dW = full([bty_zn_wMC_dW, btx_zn_wMC_dW, btt_zn_wMC_dW; ...
                          sey_zn_wMC_dW, sex_zn_wMC_dW, set_zn_wMC_dW; ...
                          pjy_zn_wMC_dW, pjx_zn_wMC_dW, pjt_zn_wMC_dW; ...
                          p0y_zn_wMC_dW, p0x_zn_wMC_dW, p0t_zn_wMC_dW]);


    results_cor_dW =full([0             , 0             , cor_MC_dW    ]);

    results_obs_dW = Nobs*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dW_j = [results_zw_dW; bar; results_zn_dW; bar; results_cor_dW; bar; results_obs_dW]; 
    clear results_zw_dW results_zn_dW results_cor_dW results_obs_dW bt* se* rj* r0* pj* p0*
    
%% Import specification

disp('running estimation for import outcomes') 

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dM, dx_MC_dM]);
    cor_MC_dM = correl_pred(2,1);
    
%Test based on naive IV
    share_wMC_nGE_dM  = share_wMC_nGE_dW(indMn==1,:);
    zn_wMC_dM = share_wMC_nGE_dM*shiftersIV;  
    
    [btt_zn_wMC_dM , set_zn_wMC_dM , rjt_zn_wMC_dM , r0t_zn_wMC_dM, pjt_zn_wMC_dM , p0t_zn_wMC_dM ] = implement_test_cluster(dytn_MC_dM, zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, critical, 0, cluster_vec);
    [bty_zn_wMC_dM , sey_zn_wMC_dM , rjy_zn_wMC_dM , r0y_zn_wMC_dM, pjy_zn_wMC_dM , p0y_zn_wMC_dM ] = implement_test_cluster(dyn_MC_dM , zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dM , sex_zn_wMC_dM , rjx_zn_wMC_dM , r0x_zn_wMC_dM, pjx_zn_wMC_dM , p0x_zn_wMC_dM ] = implement_test_cluster(dxn_MC_dM , zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dM  = share_wMC_dW_dW(indM==1,:);
    zw_wMC_dM = share_wMC_dW_dM*shiftersIV;  
    
    [btt_zw_wMC_dM , set_zw_wMC_dM , rjt_zw_wMC_dM , r0t_zw_wMC_dM, pjt_zw_wMC_dM , p0t_zw_wMC_dM ] = implement_test_cluster(dyt_MC_dM, zw_wMC_dM , share_wMC_dW_dM , shiftersIV, critical, 0, cluster_vec);
    [bty_zw_wMC_dM , sey_zw_wMC_dM , rjy_zw_wMC_dM , r0y_zw_wMC_dM, pjy_zw_wMC_dM , p0y_zw_wMC_dM ] = implement_test_cluster(dy_MC_dM , zw_wMC_dM , share_wMC_dW_dM , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dM , sex_zw_wMC_dM , rjx_zw_wMC_dM , r0x_zw_wMC_dM, pjx_zw_wMC_dM , p0x_zw_wMC_dM ] = implement_test_cluster(dx_MC_dM , zw_wMC_dM , share_wMC_dW_dM , shiftersIV, critical, 0, cluster_vec);

%Output
    
    results_zw_dM = full([bty_zw_wMC_dM, btx_zw_wMC_dM, btt_zw_wMC_dM; ...
                          sey_zw_wMC_dM, sex_zw_wMC_dM, set_zw_wMC_dM; ...
                          pjy_zw_wMC_dM, pjx_zw_wMC_dM, pjt_zw_wMC_dM; ...
                          p0y_zw_wMC_dM, p0x_zw_wMC_dM, p0t_zw_wMC_dM]);

    results_zn_dM = full([bty_zn_wMC_dM, btx_zn_wMC_dM, btt_zn_wMC_dM; ...
                          sey_zn_wMC_dM, sex_zn_wMC_dM, set_zn_wMC_dM; ...
                          pjy_zn_wMC_dM, pjx_zn_wMC_dM, pjt_zn_wMC_dM; ...
                          p0y_zn_wMC_dM, p0x_zn_wMC_dM, p0t_zn_wMC_dM]);


    results_cor_dM =full([0             , 0             , cor_MC_dM    ]);

    results_obs_dM = NM*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dM_j = [results_zw_dM; bar; results_zn_dM; bar; results_cor_dM; bar; results_obs_dM]; 
    clear results_zw_dM results_zn_dM results_cor_dM results_obs_dM share_wMC_nGE_dM share_wMC_dW_dM bt* se* rj* r0* pj* p0* 
    
 %% Revenue specification

disp('running estimation for revenue outcomes')

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dT, dx_MC_dT]);
    cor_MC_dT = correl_pred(2,1);
    
%Test based on naive IV
    share_wMC_nGE_dT  = share_wMC_nGE_dW(indTn==1,:);
    zn_wMC_dT = share_wMC_nGE_dT*shiftersIV;  
    
    [btt_zn_wMC_dT , set_zn_wMC_dT , rjt_zn_wMC_dT , r0t_zn_wMC_dT, pjt_zn_wMC_dT , p0t_zn_wMC_dT ] = implement_test_cluster(dytn_MC_dT, zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, critical, 0, cluster_vec);
    [bty_zn_wMC_dT , sey_zn_wMC_dT , rjy_zn_wMC_dT , r0y_zn_wMC_dT, pjy_zn_wMC_dT , p0y_zn_wMC_dT ] = implement_test_cluster(dyn_MC_dT , zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dT , sex_zn_wMC_dT , rjx_zn_wMC_dT , r0x_zn_wMC_dT, pjx_zn_wMC_dT , p0x_zn_wMC_dT ] = implement_test_cluster(dxn_MC_dT , zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dT  = share_wMC_dW_dW(indT==1,:);
    zw_wMC_dT = share_wMC_dW_dT*shiftersIV;  
    
    [btt_zw_wMC_dT , set_zw_wMC_dT , rjt_zw_wMC_dT , r0t_zw_wMC_dT, pjt_zw_wMC_dT , p0t_zw_wMC_dT ] = implement_test_cluster(dyt_MC_dT, zw_wMC_dT , share_wMC_dW_dT , shiftersIV, critical, 0, cluster_vec);
    [bty_zw_wMC_dT , sey_zw_wMC_dT , rjy_zw_wMC_dT , r0y_zw_wMC_dT, pjy_zw_wMC_dT , p0y_zw_wMC_dT ] = implement_test_cluster(dy_MC_dT , zw_wMC_dT , share_wMC_dW_dT , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dT , sex_zw_wMC_dT , rjx_zw_wMC_dT , r0x_zw_wMC_dT, pjx_zw_wMC_dT , p0x_zw_wMC_dT ] = implement_test_cluster(dx_MC_dT , zw_wMC_dT , share_wMC_dW_dT , shiftersIV, critical, 0, cluster_vec);

%Output
    
    results_zw_dT = full([bty_zw_wMC_dT, btx_zw_wMC_dT, btt_zw_wMC_dT; ...
                          sey_zw_wMC_dT, sex_zw_wMC_dT, set_zw_wMC_dT; ...
                          pjy_zw_wMC_dT, pjx_zw_wMC_dT, pjt_zw_wMC_dT; ...
                          p0y_zw_wMC_dT, p0x_zw_wMC_dT, p0t_zw_wMC_dT]);

    results_zn_dT = full([bty_zn_wMC_dT, btx_zn_wMC_dT, btt_zn_wMC_dT; ...
                          sey_zn_wMC_dT, sex_zn_wMC_dT, set_zn_wMC_dT; ...
                          pjy_zn_wMC_dT, pjx_zn_wMC_dT, pjt_zn_wMC_dT; ...
                          p0y_zn_wMC_dT, p0x_zn_wMC_dT, p0t_zn_wMC_dT]);


    results_cor_dT =full([0             , 0             , cor_MC_dT    ]);

    results_obs_dT = NT*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dT_j = [results_zw_dT; bar; results_zn_dT; bar; results_cor_dT; bar; results_obs_dT]; 
    clear results_zw_dT results_zn_dT results_cor_dT results_obs_dT share_wMC_nGE_dT share_wMC_dW_dT bt* se* rj* r0* pj* p0* 
    
  %% Export specification

disp('running estimation for export outcomes')

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dX, dx_MC_dX]);
    cor_MC_dX = correl_pred(2,1);
    
%Test based on naive IV
    share_wMC_nGE_dX  = share_wMC_nGE_dW(indXn==1,:);
    zn_wMC_dX = share_wMC_nGE_dX*shiftersIV;  
    
    [btt_zn_wMC_dX , set_zn_wMC_dX , rjt_zn_wMC_dX , r0t_zn_wMC_dX, pjt_zn_wMC_dX , p0t_zn_wMC_dX ] = implement_test_cluster(dytn_MC_dX, zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, critical, 0, cluster_vec);
    [bty_zn_wMC_dX , sey_zn_wMC_dX , rjy_zn_wMC_dX , r0y_zn_wMC_dX, pjy_zn_wMC_dX , p0y_zn_wMC_dX ] = implement_test_cluster(dyn_MC_dX , zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dX , sex_zn_wMC_dX , rjx_zn_wMC_dX , r0x_zn_wMC_dX, pjx_zn_wMC_dX , p0x_zn_wMC_dX ] = implement_test_cluster(dxn_MC_dX , zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dX  = share_wMC_dW_dW(indX==1,:);
    zw_wMC_dX = share_wMC_dW_dX*shiftersIV;  
    
    [btt_zw_wMC_dX , set_zw_wMC_dX , rjt_zw_wMC_dX , r0t_zw_wMC_dX, pjt_zw_wMC_dX , p0t_zw_wMC_dX ] = implement_test_cluster(dyt_MC_dX, zw_wMC_dX , share_wMC_dW_dX , shiftersIV, critical, 0, cluster_vec);
    [bty_zw_wMC_dX , sey_zw_wMC_dX , rjy_zw_wMC_dX , r0y_zw_wMC_dX, pjy_zw_wMC_dX , p0y_zw_wMC_dX ] = implement_test_cluster(dy_MC_dX , zw_wMC_dX , share_wMC_dW_dX , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dX , sex_zw_wMC_dX , rjx_zw_wMC_dX , r0x_zw_wMC_dX, pjx_zw_wMC_dX , p0x_zw_wMC_dX ] = implement_test_cluster(dx_MC_dX , zw_wMC_dX , share_wMC_dW_dX , shiftersIV, critical, 0, cluster_vec);

%Output
    
    results_zw_dX = full([bty_zw_wMC_dX, btx_zw_wMC_dX, btt_zw_wMC_dX; ...
                          sey_zw_wMC_dX, sex_zw_wMC_dX, set_zw_wMC_dX; ...
                          pjy_zw_wMC_dX, pjx_zw_wMC_dX, pjt_zw_wMC_dX; ...
                          p0y_zw_wMC_dX, p0x_zw_wMC_dX, p0t_zw_wMC_dX]);

    results_zn_dX = full([bty_zn_wMC_dX, btx_zn_wMC_dX, btt_zn_wMC_dX; ...
                          sey_zn_wMC_dX, sex_zn_wMC_dX, set_zn_wMC_dX; ...
                          pjy_zn_wMC_dX, pjx_zn_wMC_dX, pjt_zn_wMC_dX; ...
                          p0y_zn_wMC_dX, p0x_zn_wMC_dX, p0t_zn_wMC_dX]);


    results_cor_dX =full([0             , 0             , cor_MC_dX    ]);

    results_obs_dX = NX*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dX_j = [results_zw_dX; bar; results_zn_dX; bar; results_cor_dX; bar; results_obs_dX]; 
    clear results_zw_dX results_zn_dX results_cor_dX results_obs_dX share_wMC_nGE_dX share_wMC_dW_dX bt* se* rj* r0* pj* p0* 
    
%% save output

    results_dW(:,:, j) = results_dW_j;
    results_dM(:,:, j) = results_dM_j;
    results_dT(:,:, j) = results_dT_j;
    results_dX(:,:, j) = results_dX_j;

    clear share_wMC_nGE_dW share_wMC_dW_dW zw* zn* *_j

    save(save_estimation_path + "est_Tab1a.mat", 'results_dW', 'results_dM', 'results_dT', 'results_dX', 'test_cluster_grid', 'trade_target');

toc    
end

%% Visualizing results

load(save_estimation_path + "est_Tab1a.mat")

j=1;
Tab1 = results_dW(:,:, j) ;

Tab1_panelA = [Tab1(11,3), Tab1(6,3), Tab1(1,3); 
                0        ,  Tab1(7,3), Tab1(2,3);
                0        ,  Tab1(9,3), Tab1(4,3)];
